home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_Jn.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.8 KB  |  191 lines
  1. /* specfunc/bessel_Jn.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_pow_int.h>
  26. #include "bessel.h"
  27. #include "bessel_amp_phase.h"
  28. #include "bessel_olver.h"
  29. #include <gsl/gsl_sf_bessel.h>
  30.  
  31.  
  32.  
  33. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  34.  
  35.  
  36. int gsl_sf_bessel_Jn_e(int n, double x, gsl_sf_result * result)
  37. {
  38.   int sign = 1;
  39.  
  40.   if(n < 0) {
  41.     /* reduce to case n >= 0 */
  42.     n = -n;
  43.     if(GSL_IS_ODD(n)) sign = -sign;
  44.   }  
  45.  
  46.   if(x < 0.0) {
  47.     /* reduce to case x >= 0. */
  48.     x = -x;
  49.     if(GSL_IS_ODD(n)) sign = -sign;
  50.   }
  51.  
  52.   /* CHECK_POINTER(result) */
  53.  
  54.   if(n == 0) {
  55.     gsl_sf_result b0;
  56.     int stat_J0 = gsl_sf_bessel_J0_e(x, &b0);
  57.     result->val = sign * b0.val;
  58.     result->err = b0.err;
  59.     return stat_J0;
  60.   }
  61.   else if(n == 1) {
  62.     gsl_sf_result b1;
  63.     int stat_J1 = gsl_sf_bessel_J1_e(x, &b1);
  64.     result->val = sign * b1.val;
  65.     result->err = b1.err;
  66.     return stat_J1;
  67.   }
  68.   else {
  69.     if(x == 0.0) {
  70.       result->val = 0.0;
  71.       result->err = 0.0;
  72.       return GSL_SUCCESS;
  73.     }
  74.     else if(x*x < 10.0*(n+1.0)*GSL_ROOT5_DBL_EPSILON) {
  75.       gsl_sf_result b;
  76.       int status = gsl_sf_bessel_IJ_taylor_e((double)n, x, -1, 50, GSL_DBL_EPSILON, &b);
  77.       result->val  = sign * b.val;
  78.       result->err  = b.err;
  79.       result->err += GSL_DBL_EPSILON * fabs(result->val);
  80.       return status;
  81.     }
  82.     else if(GSL_ROOT3_DBL_EPSILON * x > (n*n+1.0)) {
  83.       int status = gsl_sf_bessel_Jnu_asympx_e((double)n, x, result);
  84.       result->val *= sign;
  85.       return status;
  86.     }
  87.     else if(n > 50) {
  88.       int status = gsl_sf_bessel_Jnu_asymp_Olver_e((double)n, x, result);
  89.       result->val *= sign;
  90.       return status;
  91.     }
  92.     else {
  93.       double ans;
  94.       double err;
  95.       double ratio;
  96.       double sgn;
  97.       int stat_b;
  98.       int stat_CF1 = gsl_sf_bessel_J_CF1((double)n, x, &ratio, &sgn);
  99.  
  100.       /* backward recurrence */
  101.       double Jkp1 = GSL_SQRT_DBL_MIN * ratio;
  102.       double Jk   = GSL_SQRT_DBL_MIN;
  103.       double Jkm1;
  104.       int k;
  105.  
  106.       for(k=n; k>0; k--) {
  107.     Jkm1 = 2.0*k/x * Jk - Jkp1;
  108.     Jkp1 = Jk;
  109.     Jk   = Jkm1;
  110.       }
  111.  
  112.       if(fabs(Jkp1) > fabs(Jk)) {
  113.         gsl_sf_result b1;
  114.     stat_b = gsl_sf_bessel_J1_e(x, &b1);
  115.     ans = b1.val/Jkp1 * GSL_SQRT_DBL_MIN;
  116.     err = b1.err/Jkp1 * GSL_SQRT_DBL_MIN;
  117.       }
  118.       else {
  119.         gsl_sf_result b0;
  120.     stat_b = gsl_sf_bessel_J0_e(x, &b0);
  121.     ans = b0.val/Jk * GSL_SQRT_DBL_MIN;
  122.     err = b0.err/Jk * GSL_SQRT_DBL_MIN;
  123.       }
  124.  
  125.       result->val = sign * ans;
  126.       result->err = fabs(err);
  127.       return GSL_ERROR_SELECT_2(stat_CF1, stat_b);
  128.     }
  129.   }
  130. }
  131.  
  132.  
  133. int
  134. gsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double * result_array)
  135. {
  136.   /* CHECK_POINTER(result_array) */
  137.  
  138.   if(nmin < 0 || nmax < nmin) {
  139.     int n;
  140.     for(n=nmax; n>=nmin; n--) {
  141.       result_array[n-nmin] = 0.0;
  142.     }
  143.     GSL_ERROR ("domain error", GSL_EDOM);
  144.   }
  145.   else if(x == 0.0) {
  146.     int n;
  147.     for(n=nmax; n>=nmin; n--) {
  148.       result_array[n-nmin] = 0.0;
  149.     }
  150.     if(nmin == 0) result_array[0] = 1.0;
  151.     return GSL_SUCCESS;
  152.   }
  153.   else {
  154.     gsl_sf_result r_Jnp1;
  155.     gsl_sf_result r_Jn;
  156.     int stat_np1 = gsl_sf_bessel_Jn_e(nmax+1, x, &r_Jnp1);
  157.     int stat_n   = gsl_sf_bessel_Jn_e(nmax,   x, &r_Jn);
  158.     int stat = GSL_ERROR_SELECT_2(stat_np1, stat_n);
  159.  
  160.     double Jnp1 = r_Jnp1.val;
  161.     double Jn   = r_Jn.val;
  162.     double Jnm1;
  163.     int n;
  164.  
  165.     if(stat == GSL_SUCCESS) {
  166.       for(n=nmax; n>=nmin; n--) {
  167.         result_array[n-nmin] = Jn;
  168.         Jnm1 = -Jnp1 + 2.0*n/x * Jn;
  169.         Jnp1 = Jn;
  170.         Jn   = Jnm1;
  171.       }
  172.     }
  173.     else {
  174.       for(n=nmax; n>=nmin; n--) {
  175.         result_array[n-nmin] = 0.0;
  176.       }
  177.     }
  178.  
  179.     return stat;
  180.   }
  181. }
  182.  
  183. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  184.  
  185. #include "eval.h"
  186.  
  187. double gsl_sf_bessel_Jn(const int n, const double x)
  188. {
  189.   EVAL_RESULT(gsl_sf_bessel_Jn_e(n, x, &result));
  190. }
  191.